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A successive continuation method for locating connecting orbits in parametrized 
systems of autonomous ODEs was considered in [9]. In this paper we present an 
improved algorithm for locating and continuing connecting orbits, which includes a 
new algorithm for the continuation of invariant subspaces. The latter algorithm is of 
independent interest, and can be used in different contexts than the present one. 
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1 Introduction. 

Homoclinic and heteroclinic orbits, also called connecting orbits, are trajectories connecting 
equilibrium points of a system of autonomous ordinary differential equations. Computation 
of connecting orbits is becoming increasingly important, both in dynamical systems research, 
such as understanding chaotic dynamics, and in a variety of applied problems, including wave 
propagation in combustion models, chemical reactions, neuronal interactions, solitary waves 
in fluid, solitons in nonlinear optical fiber, and communication processes in living cells, to 
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name a few. The corresponding numerical problem is that of finding solutions (u(t),X) of 
the system of autonomous ODEs 



u'(t) - f(u(t),X) = 0, «(•), /(•, •) e W ,X e W 

lim u(t) = uq, lim u(t) = u\. 



t— ¥ — OO t— ¥+00 



(1) 

(2) 



Most algorithms for the numerical analysis of connecting orbits reduce (1), (2) to a 
boundary value problem on a finite interval using linear or higher order approximations of 
stable and unstable manifolds near uq and u\, respectively. See recent papers by Champneys, 
Kuznetsov, Sandstede [4], by Doedel, Friedman, Kunin [9] and Moore [13] for the history of 
the question and the bibliography. Note that in the last work an alternative approach was 
used based on the arclength parametrization, instead of using time t as a parameter. 

The algorithms in [4] use a version of Beyn's continuation algorithm based on projection 
boundary conditions ([1], [2]). They were implemented in a set of routines, HomCont, which 
are currently part of AUT097 [7]. HomCont has capabilities for detailed bifurcation analysis 
of homoclinic orbits and some bifurcations of heteroclinic orbits. It has limited capabilities 
for locating connecting orbits, namely, a simplified version of the algorithm in [9]. 

The algorithms in [9] have their primary focus on locating connecting orbits and use a 
modification of a continuation algorithm based on projection boundary conditions (Friedman, 
Doedel [11]). They were implemented in an experimental code based on AUT094 [8]. 

In order to have a well posed problem, it is necessary for the boundary conditions to be 
sufficiently smooth with respect to parameters. Both in [4] and [9], the boundary conditions 
are defined with respect to bases of stable or unstable eigenspaces of f u (uo, A) an d f u ( u i, 
The approach in [4] is to compute an orthonormal basis in the appropriate eigenspace, at 
each pseudo arclength continuation step, using a "black box" routine based on the real Schur 
factorization and then to adapt this basis to be smooth with respect to parameters, using 
a technique due to Beyn [2, App. C] which amounts to the solution of a linear system of 
the dimension of the eigenspace in question. The approach in [9] is to compute initially an 
orthonormal basis in the appropriate eigenspace via the real Schur factorization and then 
to continue the real Schur factorization equations (as a part of boundary conditions). At 
the same time precise convergence of the algorithm in [9] is not clear, and it is somewhat 
cumbersome to use. In some recent work, [6], Dieci and Eirola provide a general differential 
equations framework for continuation of the block Schur factorization as well as other matrix 
factorizations. Reference [6] includes a comprehensive set of references for smooth matrix 
factorization for parameter dependent matrices. 

In this paper we present an improved algorithm for locating and continuing connecting 
orbits, which includes a new algorithm for the continuation of invariant subspaces (CIS). 
This CIS algorithm is based on iterative refinement techniques originally due to Stewart 
[14], and later revisited by Demmel [5]. We provide some new twists to these techniques: (i) 
we justify these iterative refinement techniques using the differential equations which model 
continuation of block Schur forms, and (ii) we make use of these differential equations to 
obtain an accurate approximation of the relevant invariant subspace. 

In the end, the new algorithm is more efficient than the algorithms in [4] and [9] and 
is very robust. In particular, it provides several possible safeguards against fast variation 
of eigenvalues. It has been implemented in an experimental code based on AUT097 which 
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is essentially a modification of the HomCont part of AUT097 to include the algorithm in [9] 
for locating and continuing connecting orbits and the CIS algorithm, while preserving the 
bifurcation analysis part of HomCont. 



2 An improved algorithm for locating and continuing 
connecting orbits. 

Assume, for simplicity of notation, that the fixed points u and U\ are hyperbolic, and the 
eigenvalues of f u (u , X) and f u (ui,X), respectively, satisfy 

Re/V < ••• < Re ^o,n +i < < A*o,i < Re ^o, 2 < ••• < Re^o,r.o' 
Refi 1;1 < ... < Re/i lni < < Re/i lni+1 < ... < Re/i ln . 

In this paper, we will assume that the matrices f u (uo,i, A) are smooth functions of A (for A 
in an appropriate subset of M" A . 

The method extends to the case [i 1 = 0, as in [4]. It also extends to the cases of complex 
and multiple [i 1 by a simple modification of Step 0, eq. (11) below, of the algorithm 
(see [9, Section 4.3] for a computational example). The algorithm requires evaluation of 
various projections associated with the eigenspaces of f u (u , A) and f u {ui, A). Initially these 
projections are constructed using the real Schur factorizations [12] 

f u (u , A) = Q T Qq, f u {ui, A) = QiTiQ\. 

The first factorization has been chosen so that the first no columns qo,i, ■■■,Qo,n of Qo 
form an orthonormal basis of the right invariant subspace So of f u { u o,X), corresponding 
to fj, 01 , fj, , and the remaining n — riQ columns go,n +i i ■ ■ ■ ■> Qo,n °f Qo form an orthonormal 
basis of the orthogonal complement Sq. Similarly, the first ri\ columns q^i, q\^ n of Qi 
form an orthonormal basis of the right invariant subspace Si of f u (ui,X), corresponding to 
■■■i Ll i,n 1 i an( ^ the remaining n — ni columns ?i ini +i, ■■■,q\, n °f Qi form an orthonormal 
basis of the orthogonal complement S^. In the algorithm below the matrices Qq(X) and 
Qi(A) are assumed to be computed at each continuation step by a "black box" routine, 
described in Section 3, which ensures their continuity. 

The approximate finite interval problem is to find a branch of solutions (u(t), A, u , Ui, T), 
u £ ^([0, 1],K"), A G R" A , provided n A = n - (n + n x ) + 2; n A > 0, u ,Ui e K", T > 0, 
is the length of the time interval, for some small eo, e i > 0, of the time-scaled differential 
equation 

u'(t) - Tf(u(t),X) = 0, 0<t<l, (3) 
subject to left boundary conditions 

(u(0) - u ) ■ q ,n +i(u ,X) = 0, % = 1, ...,n - n , (4) 
\u(0)-u \ = e , (5) 

right boundary conditions 

(tt(l) -ui) ■ qi :ni+i (ui, X) = 0, % = 1, ..., = n - n x . (6) 
|m(1)-Mi| = ei (7) 
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stationary state conditions 

/K,A) = 0, (8) 
f(u u X) = 0. (9) 

Remark 1 Initially, we perform time integration to obtain a, typically, crude orbit with 
initial point u(0) G So but the terminal point u(l) ^ Si, in general. Hence defined by 

Ti = (u(l) - ui) ■ gi, ni +i(«i, A)/ei, i = 1, ...,n T = n - m, (10) 

are, in general, nonzero, and the initial connecting orbit on the branch of connecting orbits 
is found via a sequence of homotopies that locate successive zero intercepts of the Tj in (10). 
In each homotopy step we compute a branch, i.e., a one- dimensional manifold, of solutions. 
For this we must have n c — n v = n — 1, where n c is the number of constraints, and n v is the 
number of scalar variables. We keep u(l) free, while u (0) is allowed to vary on the surface 
of the sphere in So of radius : (i) according to equation (12) below at steps through no, 
when X is fixed and q in (12) play the role of the control parameters; and (ii) according to 
equation (5) at steps n + 1 through n + n\, when X varies. 

Let So,*:, k = 1, no, be the right invariant subspace of f(uo, Ao) corresponding to the 
eigenvalues /i l5 ji Q k . Then the first k columns qo,i, ■■■,qo,k of Qq form an orthonormal 
basis of So,* and the remaining n — k. columns qo^+i, Qo,n of Qo form an orthonormal basis 
of the orthogonal complement S^. 

1. Initialization. 

Step 0. Initialize the problem parameter vector A, and set the algorithm parameters e 
and T to small, positive values, so that u(t) is approximately constant on [0, T]. Set 

u(t) = u + e ci?oi, < t < 1, (11) 

or 

u(t) = u + eoci^oie 116 ^ ' 1 *, < * < 1, Re u 0)1 > 0, 
ei = \u(l) — u\\ , c\ = 1, or —1 and c^ = ... = c no = 0. 

2. Locating a connecting orbit, A is fixed. 

Step 1. Time integration to get an initial orbit. Compute a solution branch to the 
system (3), (4), (10), (7), and 

{u(0)-u Q )-q Qi i(uQ,X)/eo-Ci = 0, i = 1, n c = n , (12) 

in the direction of increasing T, until u(l) reaches an ei-neighborhood of U\, for some 
ei > 0. Scalar variables are T, e\ e M, r e M" - " 1 . There are n differential equations 
with n c = 2n — ri\ + 1 constraints and n„ = n — ni + 2 scalar variables, and hence 
n c — n v = n — 1. In practice one typically continues until ei stops decreasing, its value 
being not necessarily small. 

Step k, k = 2,...,n (for n > 1 ). Compute a branch of solutions to the system 
(3)-(5), (10), (7), (12) to locate a zero of, say, t^-i (while Ti,...,Tfc_ 2 = 0, fixed). 
Free scalar variables are ei,Ci, ...,Ck,Tk-i, ...,r„_ ni . There are n differential equations 
with n c = 2n — ri\ + 2 constraints and n„ = n — ri\ + 3 scalar variables, and hence 
n c — n v = n — 1. 
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3. Locating a connecting orbit, A varies. 

Step k, k = no + 1, ...,no + n\ = n — n\ + 1. Compute a branch of solutions to the 
system (3)-(5), (10), (7)-(9) to locate a zero of, say, Tfc_i(while ri, ...,Tk-2 = 0, fixed). 
Free scalar variables are ei, r^-i, r„_ ni , Ai, Xk- no G M, iio, ?/i G M n . There are n 
differential equations with n c = 4n — n — n\ + 2 constraints and n„ = 3n — n — ri\ + 3 
scalar variables, and hence n c — n v = n — 1. 

4. Increasing the accuracy of the connecting orbit. 

Compute a branch of solutions to the system (3)- (9) in the direction of decreasing ei 
until it is 'small'. Free scalar variables are ei, T, Ai, A„ A _i 6 R, n ,Mi G R™. As 
before, n c = 4n — no — ni + 2, n„ = 3n — no — ni + 3. 

5. Continue the connecting orbit. 

Compute a branch of solutions to the system (3)- (9). Free variables are the (real) scalar 
T, and Ai, . . . , A„ A , and the vectors u , u\ G R n . As before, n c = 4n — n — n\ + 2, 
n„ = 3n — n — ni + 3. Alternatively, a phase condition 



may be added if T is kept fixed and e and ei are allowed to vary. Here q(t) is a 
previously computed orbit on the branch. 

Remark 2 In principle, our algorithm of continuation of invariant subspaces of f u (uo, X) 
and f u (ui,X) in Section 3 breaks down only if two eigenvalues, one associated with the sub- 
space being continued and one not, approach the same point on the imaginary axis (one 
from the left and one from the right). In this case the algorithm should stop anyway since a 
bifurcation is being approached. 

3 Continuation of Invariant Subspaces. 

Let A(X) G M" x " denote one of the following: f u (uo,X), f u (ui,X). The basic continuation 
algorithm requires at each pseudo arclength continuation step computation of a right in- 
variant (typically, stable or unstable) m-dimensional subspace S(X) of A(X). In general, the 
function A is smooth in A (say, differentiable) , and it is important that 5(A) be also smooth, 
as otherwise convergence difficulties can be expected. 

In this section we show how to constructively obtain smooth and orthogonal Q(X) = 



[Qi(A)|Q 2 (A)] G R nxn , Qi(A) e M nxm , Q 2 (X) G R«x(«-™), g0 that gpan and 



(32(A) span the orthogonal complement S(X)- 1 -. Moreover, if we let Sk(X), k = l,...,m, be 
the right invariant subspaces of A(X) corresponding to the first k eigenvalues u^...,^ of 
A, then the first k columns qi, ...,qk of Q(X) form an orthonormal basis of Sk(X), and the 
remaining n — k columns form an orthonormal basis of the orthogonal complement Sk{X) L . 

To justify our construction, we should recall that in our continuation procedure we 
parametrize a solution branch in terms of so called pseudo arclength; let s denote the pseudo 




(13) 
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arclength variable 1 . Thus, both fixed points u and ui, see (8)-(9), as well as the param- 
eter(s) A are smooth functions of s. The matrix valued function A : A G M nA — >■ M" x " can 
thus be viewed as a smooth function from s e R — > R" x ™. As a consequence, we can and 
will think of invariant subspaces' continuation with respect to the scalar pseudo arclength 
variable s. For this reason, we will abuse of notation and freely write A(s) for A(X). In what 
follows, a "' " will indicate differentiation with respect to s. 

Remark 3 If n\ = 1 in (1), then it is conceivable that one may be able to perform contin- 
uation with respect to A itself, rather than reparametrizing the problem by arclength. 

The basic issue is the following: "Suppose that initially we have the (real) block Schur 
factorization 

A(0) = Q(0)T(0)Q T (0), Q(0) = [Qi(0)|Q 2 (0)], (14) 
where T(0) is block upper triangular (Tu(0), i = 1,2, are not required to be triangular) 



T(0) 



Tn(o) r 12 (o) 
o r 22 (o) 



the columns of Qi(0) span an invariant subspace 5(0) of A(0), and the columns of Q 2 (0) 
span the orthogonal complement 5(0) ± . We want to obtain a block Schur factorization for 
the matrix A(s), close to A(0), exploiting (if possible) the work already done to obtain the 
block Schur form for A(0)" . This problem fits within the general framework developed in [6]. 
Suppose that the matrix A(s) has two groups of eigenvalues, Ai(s) and A 2 (s), which stay 
disjoint for all s around 0. Then, in an interval about s = there is the smooth factorization 

A(s) = Q(s)T(s)Q T (s), Q(s) = [Qi( S )|Q 2 (s)], (15) 



where T(s) is in block Schur form T(s) 



. Here, Tn has eigenvalues Ai 



T n (s) T 12 (s) 
T 22 (s) 

and T 22 has eigenvalues A 2 . A constructive procedure to obtain the factorization QTQ T of 
A can be based upon differential equations models. As in [6], we differentiate the relations 
A = QTQ T and Q T (s)Q(s) = /, let H := Q T Q, and obtain 

f = Q T AQ + TH — HT , (16) 
Q = QH. (17) 

Now we use triangularity of T and the fact that H must be skew symmetric, H T = —H; 
we partition H in the same way as T, and then can determine H\2 by solving the Sylvester 
equation 

T 22 H* - H*T n = (Q T AQ) 12 . (18) 

The blocks H n and i? 22 are not uniquely determined, and we may set them to (in any case, 
they must be chosen skew symmetric). Thus, in principle, we can solve (16)-(17) subject to 



1 Typically, this is the general procedure of most practical continuation algorithms; e.g., this is the strategy 
implemented in AUTO. 
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initial conditions (ICs) obtained from the factorization of A(0), in order to obtain a smooth 
path of block Schur factorizations. 

We can accumulate the transformations in such a way that we are always looking for 
corrections close to the identity. To be more precise, we can rewrite 



Q(s) = Q(0)U(s) , with [7(0) 



(19) 



and use this in (17), thereby obtaining a differential equation for U, U = UH (notice that 
H is the same as before). We now look for exact solutions to the U differential equation. 

U n {s) U 12 (s) 
U 21 (s) U 22 (s) 

Since U(0) = I, there is an open interval about where we can require that U\ has the 

U n - Next, we define 



Partition U(s) = [Ui(s), U 2 (s) 
Since U(0) = 
structure U\ 



with same block dimensions as T. 



U 21 U^ 



Y(s) := U 21 { 8 )UJ{8), 



and use the orthogonality relation UjU\ = I to obtain U\ 



I 

Y 



(20) 

(/ + Y T Y)- l l 2 . In a 



similar way, for U 2 (s) we use U 2 U 2 = I and UfU 2 = 0, to eventually obtain 

U(s) = [( / (s) ) (/ + Y T Y)~ 1 ^ 2 , ( ) (I + YY T )-^ . 



(21) 



Thus, we need to find the matrix Y G R( n - m )* m in (21). Define E(s) by 



Q T (0)A(s)Q(0) = Q T (0)[A(0) + (A(s) - A{0))]Q{0) =: T(0) + E(s) 



Tn T12 
E 2 \ T 22 



■ (22) 



Now we substitute Q(s) given by (19), (21) and A(s) obtained from (22) into the invariant 
subspace relation: 

QUs)A(s)Q 1 (s) = 0, (23) 



to obtain the following algebraic Riccati equation for Y : 

f 22 Y - YT n = -E 2l + YT 12 Y, 



(24) 



or 

F(Y) = , F(Y) := f 22 Y - YT n + E 2l - YT 12 Y. (25) 

Remark 4 If we think 0/5(0) and Q(0) as approximations to S(s) and Q(s), then we may 
interpret the form (21) as well as the resulting Riccati equation (24) as an iterative refinement 
technique to improve the accuracy of computed invariant subspaces. This is the viewpoint in 
the works of Stewart [14], Dongarra, Moler and Wilkinson [10], Chatelin [3], and Demmel 
[5]. However, we prefer to think of (21) as the exact solution of the differential equation 
U = UH (see (17) and (19)). We will exploit this fact later. 
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How to solve (24). The following two iterative methods have been often advocated to 
solve the Riccati equation (24) (or (25)): 



1. The iteration [14]: 

f 22 A k - A k f u = , Y k = A k + Y k _ x , (26) 

with Y Q = 0, k = 1,2, . . . . 

2. The Newton iteration [5]: 

(f 22 - n_if 12 )A fc - A k (f u + f 12 y fc _!) = , Y k = A k + y fc _! , (27) 

with y = 0, A; = 1,2, . . . . 

Therefore, we need to solve a Sylvester equation in the inner loop of the iterative refine- 
ment. This can be effectively done by using LAPACK routines. As shown in the convergence 
analysis for the iterations (26) and (27) by Stewart [14] and Demmel [5], respectively, if we 
let 

k = — - — , (28) 

sep 2 (T n ,T 22 ) 

then under the assumptions k < 1/4 and k < 1/12, the iterations (26) and (27) converge, 
linearly and quadratically, respectively. In (28), || • is the Frobenius norm of a matrix. 

The parameter n can be interpreted as follows. Its numerator, HT^HfII-EW ||f measures the 
quality of the initial approximate invariant subspace: it will be small when the approximation 
is good, and the factor H-E^iHf will be zero if and only if the initial approximation is in fact 
correct. The function sep(T 11 ,T 22 ) in the denominator is the smallest singular value of the 
operator which maps Y to T 22 Y — YT n , and it measures the separation of the spectra of 
Tn and T 22 . If sep(T 11 ,T 22 ) is small, it means that some eigenvalues of T u and T 22 can be 
made to merge with small changes in Tu; this means that the invariant subspaces belonging 
to the two parts of the spectrum are unstable and hard to compute. Thus n will be small if 
we start with a good initial approximate invariant subspace and if the eigenvalues associated 
with that subspace are well separated from the remainder of the spectrum. Oversimplifying, 
both algorithms converge if (i) the spectra of T n (0) and T 22 (0) are far enough apart, and 
(ii) the perturbation E(s) of T(0) is small enough. 

Remark 5 We point to some additional safeguards, which ensure proper performance of the 
algorithms (26) and (27). 

1. Note that from (22) we have \\E 21 \\ F < \\E\\ F < \\A(s) - A(0)\\ F . Hence, from (28), if 

a\\A(s)-A(0)\\F<- r a= 5"^. , (29) 
4 sep 2 (T n ,T 2 2) 

then k < 1/4. Therefore we are guaranteed that for any matrix X G W ixn with 
a\\X — A(0)\\p < \ we can find the invariant subspaces S and S 1 - by Algorithm (26), 
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say. And no eigenvalue of A \s can 'merge' with an eigenvalue from A \ s ±, where A \s 
denotes the restriction of A to S, etc.. In other words, the eigenvalues of A \$ and 
A \ s ± remain separated for all matrices in a ball around A(0) of radius l/(4a). Thus, 
employing the safeguard (29) ensures that the iterations (26) and (27) will diverge only 
when a small perturbation A(s) of A(0) will make an eigenvalue from S(s) and an 
eigenvalue from S(s) 1 - coalesce. 

2. The quantity H-E^iHf / 'sep(Tn, T22) can be interpreted as the tangent of the angle between 
the subspaces spanned by Qi(0) and Q\(s). Hence the convergence of our algorithms 
implies that this angle always stays between and n/2. If one wants to monitor this 
angle, a convenient measure is the sine of the angle given [12] by 

dist^o), q^)) = ^i-o-i iQ {omo^sY). 

Substituting into this equation the formula for Qi(s) gives 



dist^o), ^(s)) = vT^mr 



1/2 _ ||^|| 2 



1 + IIF" 2 



(30) 
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The previous discussion has been motivated by the standard linear algebra viewpoint, 
that is of "how to refine the initial trivial estimate U = P. In particular, this gave us the 
initial conditions Y = for the iterations (26) and (27). However, from the point of view 
of continuation, this is usually not the best strategy, since it amounts to starting with the 
old solution as initial guess. We can do better by using the differential equation formulation 
and taking an Euler approximation to the solution. 

Recall that the matrix Y(s), solution of (24), is related to [/by (20): Y(s) = U 2 i(s)U{ 1 1 (s). 

11 1 f [U\, U 2 ] solves the differential equation U = UH (see (17)- 



Here, U\ - 
(19)). That is 



U, 



21 



and U 



U l = U l H ll + U 2 H 2l , U 2 
Recall that the factor H 2 i(s) is determined by (18 



= -KH^ + U 2 H 22 . (31) 

but there is freedom insofar as the choice 
of Hn and H 22 . We now show that Y(s) is unaffected by the choices made for Hn (and ^22)- 
In particular, this implies that the range of s-values guaranteeing invertibility of Un(s), and 
hence the representation (21), is unaffected by the choice of Hn and H 22 in (31), and hence 
there is little reason not to set H u and H 22 both equal to when we represent U as in (21). 

Lemma 1 Let U(s) = [Ui(s), U 2 (s)] be the solution of (31) with U(0) = I, obtained by 
setting Hn(s) = and H 22 (s) = in (31), while determining H 2 \(s) by (18). Let U(s) = 
[Ui, U 2 ] be the solution of (31) with U(0) = I, and Hn and H 22 nonzero skew symmetric 
matrices. Then, we have 

U^s) = C/i(s)C(s) , 



where C(s) E R mXTO is orthogonal. Therefore, partitioning U\ 
U u we have Y(s) = U^U^is) = U 2 i{s)U^ (s) . 



U, 



21 



and similarly for 
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Proof. Let p(X, s) be the characteristic polynomial of the matrix A(s). Then, we have the 
factorization p(X, s) = p\(X, s)p 2 (X, s) where the factors p\ and p 2 have no common root, and 
the roots of p\ give A 1; while the roots of p 2 give A 2 . Thus, we have 

p 1 (Q T (0)A(s)Q(0), 3)^(3) = 0, and 

Pl (Q T (0)A(s)Q(0), 3)0,(3) = . 

Therefore, since iank(pi(A(s) , sj) = n — m, then Ui(s) = Ui(s)C(s), where C(s) is orthog- 
onal (since tiff (s)Ui(s) = I). ■ 

Now, because of Lemma 1 we can just integrate the differential equations (31) with Hn 
and H 22 equal to 0: 

[U u U 2 ] = [U 2 H 21 , , U(0) = 1. 

So, the idea is to approximate the solution of these differential equations from to s in order 
to get an approximation V\ to U\(s) and thus obtain an approximation to Y(s) from V^iV^ 1 , 

V n 



where we have partitioned V\ 



21 



. We use a forward Euler step to obtain V\. 



1. Initialization. 

Set U(0) = I and obtain H 2 \(0) from solving the Sylvester equation 

T 22 (0)# 21 (0) - # 21 (0)T n (0) = - (0, /) Q(0) T i(0)Q(0) (/, 0) T . (32) 

2. Euler step. 

Let V\ = Ui(0) + sU 2 (0)H 2 i(0). Obtain initial guess for the Riccati equation (24) (or 
(25)): 

Y = V 21 V n \ (33) 

Remark 6 We now look at how this new initial guess Yq in (33) impacts convergence of the 
iterations (26) and (27). We also make some comments on expense. 

1. Observe that the value of Yq in (33) is nothing but Yq = sH 2 i(0). It is easy to verify 
that this is precisely the same approximation we would have obtained by using a forward 
Euler step to approximate the solution at s of the differential equation satisfied by Y . 

2. By our construction, we see that Y in (33) is a second order approximation (in s) to 
the exact solution Y(s). More precisely, from Taylor expansion at s = of the exact 
solution Y(s), we immediately get that Y(s) = sY(0) + 0(s 2 ) = sH 21 (0) + 0(s 2 ). On 
the other hand, the initial guess Yq = is only an O(s) approximation to Y(s). 

3. In practice, in (32), we do not have a close expression for ^4(0), but we can replace it 
by the difference quotient (l/s)(A(s) — A(0)). With this choice, the value ofY (defined 
as in (33)) turns out to be the solution of the Sylvester equation 

T 22 (0)Y -Y T n (0) = -E 21 . (34) 

Since A(0) = (l/s)(A(s) - A(0)) + 0(s 2 ), the value of Y in (34) is still a 0(s 2 ) 
approximation to the exact Y(s). 



2 Use of the Euler step is an accepted standard in continuation algorithms, and it is the strategy imple- 
mented in AUTO. 
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4- The estimates (and proof ) of convergence for (26) and (27) relied on the value of n in 
(28) to be sufficiently small. Much the same arguments can be used for the refined guess 
in (33), by the following simple modification. In U(s) in (21), we let Y(s) = Y + AY 
and use this form ofU in the invariant subspace relation (23). So doing, we obtain the 
new Riccati equation 

(f 22 - Y f 12 ) AY - AY(f n + f 12 F ) = -F{Y ) + AYT l2 AY. (35) 

Now, it is a simple verification that iterative solution of (35) by either of the two itera- 
tions (26) or (27), appropriately reformulated for the unknown AY, produces precisely 
the same sequence as the original iterations on the unknown Y had they been started 
with initial guess Y in (33). As a consequence, the convergence results we quoted after 
(28) hold unchanged, except that we have a new k value, call it k : 

« = l|ri2 ^| F( ^ 0) l lF , f n = f u + f 12 Y , f 22 = f 22 - Y f 12 . (36) 
sep 2 (T n ,T 22 ) 

We now look at the ratio 

k = \\F(Y )\\ F f sep(f u ,f 22 ) 
« ||^2i \\f ysep(T u ,T 22 ) 

By our previous remark, F(Y ) is close to at second order in s, whereas E 2X is 
only first order close to 0. To compare sep(T u ,T 22 ) with sep(T u ,T 22 ), we realize that 
sep(Tn,T 22 ) is a second order approximation to sep(Tn(s) ,T 22 (s)) , whereas sep(Tn,T 22 ) 
is only a first order approximation to the same quantity. Therefore, in first approxi- 
mation, ( | f^ n '^ 22 | ) is 1 + 0(s / sep(T u (s) ,T 22 (s))) . To summarize, we have that 
V sep(i ii ,i22) J 

k/K = 0(s), (37) 

with an obvious improvement in the radius of convergence for both iterations (26) and 
(27). 

5. By creating the initial guess Yq from (33) we need to solve the Sylvester equation (32), 
or (34)- Typically, this involves Schur reduction of the matrices Tn(0) and T 22 (Q), 
which is somewhat expensive. However, suppose we use the Newton iteration (27); 
then, at each iteration step we need to Schur factor the matrices T u + F fc _ 1 T 12 and 
T 22 — Ti 2 Y. At convergence of Newton's method, these are precisely the matrices we 
will need for the next continuation step in (32). Therefore, no extra factorizations are 
needed in this case. For a practical comparison of the iterations for the two initial 
guesses (33) and Yq = 0, we refer to the next section. 

Remark 7 An alternative to our approach for continuing orthogonal invariant subspaces 
could be easily obtained as follows. Let Q(0) : Q T (0)A(0)Q(0) be a Schur factorization 
of A(0), Q(0) = [Qi(0), $2(0)] as usual, so that (say) Qi(0) spans the invariant subspace 
relative to the eigenvalues with positive real parts. Next, consider A(s) which has the same 
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number as A(0) of eigenvalues with positive real part. Then, one may think to take the 
ordered Schur factorization of the matrix A(s) : P T A(s)P, and partition P = [Pi, P2] so 
that Pi spans the invariant subspace relative to the eigenvalues with positive real part of A(s). 
In general, Pi is not smooth (i.e., it is not true that Pi = Qi(0) + sQi(0) + ...). To enforce 
smoothness, we may solve an orthogonal Procrustes problem and replace Pi with PiV, where 
V is the orthogonal polar factor o/P 1 T Qi(0) (see [12]). Essentially, this is the approach used 
by Beyn, [2], and then implemented in HomCont, [4]. However, we believe that the approach 
we have adopted is preferable to the one we just outlined. For one thing, the latter is usually 
more expensive than our approach. Moreover, unlike the one we just outlined, our approach to 
continuation of invariant subspaces is more in tune with the original continuation problem: 
continuation of the subspaces influences the continuation step, and using first derivative 
information (as we did to obtain (33)) is bound to reflect genuine difficulties of the original 
differential equation (such as nearing a bifurcation) into the continuation algorithm. 

4 Example: Heteroclinic orbits in a 4-D singular per- 
turbation problem. 

Consider the problem of finding traveling wave front solutions to the FitzHugh-Nagumo 
equations with two diffusive variables 

v t = v xx + v(v - a)(l - v) - w, w t = 5w xx + e(v - jw), (38) 

for small positive e, while S ranging between a small and large value. For S small this is a 
singularly perturbed reaction-diffusion system. In moving coordinates, V\ = v(z), v 2 = v (z), 
Wi = w(z), w 2 = w'(z) with z = t + cx, the reduced ODE is 

v[ = v 2 , 

v' 2 = cv 2 — vi(l — vi)(vi — a) + w, (39) 

w[ = w 2 , 

w' 2 = [cw 2 — e(vi — ^Wi)]/S. 

In [9] we located a heteroclinic orbit of (39). Here we first reproduced this result with 
our new code and then located a heteroclinic orbit with S = e = 0.001, 7 = 13.23529, 
a = 0.3, and c = 0.2571271. In this case n = n\ = 2, where the relevant eigenvalues are 
H 0t i = 0.6958, yu 2 = 257.2, //^ = -0.4247, and \i X 2 = -0.06553. We then performed a two 
parameter continuation in (5, c) in the direction of increasing S, see Fig. 1 for the bifurcation 
diagram and Fig. 2 for some typical solutions in (vi,t) coordinates. 

We have implemented 4 methods to solve the Riccati equation: (i) the "Simple Iteration" 
(26) with zero initial guess, (ii) Newton's method (27) with zero initial guess, (iii) Simple 
Iteration (26) with Euler initial guess, and (iv) Newton's method with Euler initial guess. 
The numerical results agree with our theoretical results. Namely, the choice of the Euler 
initial guess is a big improvement with respect to the simpler zero guess. Specifically, for the 
same pseudo-arclength continuation step (it took 433 such steps to compute the branch), 
(iii) and (iv) required much fewer iterations to converge than (i) and (ii). Some insight in 
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Figure 1. Bifurcation diagram in (S,c) coordinates. At label 10 (<5, c) = (0.3198,0.2376), 
£/ie eigenvalues are: fi 1 = fi 2 = 0.7406, /i x : = —0.4357, and /x 1 2 = —0.06502. Za&eZ 
11 (5,c) = (0.4462,0.2265), eigenvalues are: fi Q>1 = // 02 = o'.6204, // M = -0.4409, 
and Hi2 = —0.06565. The eigenvalues are complex between the labels 10 and 11. At label 12 
(5,c) = (0.5085,0.2202), the eigenvalues are: fi 01 = 0.5098, ^ = 0.6538, ^ = -0.4436, 
and /j, 12 = —0.06612, where + \t l i,2\ = A*o,i- /a&e/ 13 (5, c) = (1.383,0), the eigen- 

values are: fj, Q 1 = —\X\i = 0.1099, fj, Q2 = — A*i i = 0.5454. And the part of the branch below 
the S-axis is symmetric to the one above the 8-axis. 

the advantages gained by use of the Euler guess is obtained by considering the typical con- 
vergence behavior of the four methods (i)-(iv): on average, method (i) required 7 iterations 
for convergence, method (ii) needed 5, method (iii) needed slightly more than 4, and method 
(iv) required less than 3 iterations for convergence. In one exceptional case, methods (i) and 

(ii) required as many as 25 iterations. Finally, (i) and (ii) failed to converge in some cases (at 
the end of the continuation), because the continuation step was not small enough, whereas 

(iii) and (iv) never failed to converge. 

Labels 10, 11, 12, 14, and 15 mark local bifurcations (of the eigenvalues), while label 
13 marks a global bifurcation (the detailed analysis of these bifurcations will be given else- 
where). For comparison, we repeated the computation of the above branch using HomCont 
(in AUT097), see Fig. 3. Note that HomCont could not continue the original branch beyond 
the global bifurcation point and instead switched the branches. We also note that by varying 
the continuation step size in our code we could select the desired branch at Label 13, while 
it was not possible to achieve this with HomCont. This confirms our theoretical insight that 
our numerical method is more in tune with the original continuation problem. 
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Figure 2. Some typical solutions in (vi,t) coordinates. 
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